############################FIRST PART################################
###################Function to get average results####################
######################################################################
Figure4f <- function(dirpath,testmethod,sim.N,numobs,TOPNMAX=300){
  setwd(dirpath)
  DASEV.TPR <-c(); TLK.TPR <-c();AFT.TPR <- c();  TWOW.TPR <-c();TWOT.TPR <-c()
  
  #cutoff=0.01
  a1_all <- c();  b1_all <- c()#dasev
  c1_all <- c();  d1_all <- c()#TLK
  e1_all <- c();  f1_all <- c()#AFT
  i1_all <- c();  j1_all <- c()#two T
  
  #cutoff=0.05
  a2_all <- c();  b2_all <- c()
  c2_all <- c();  d2_all <- c()
  e2_all <- c();  f2_all <- c()
  i2_all <- c();  j2_all <- c()
  
  #cutoff=0.1
  a3_all <- c();  b3_all <- c()
  c3_all <- c();  d3_all <- c()
  e3_all <- c();  f3_all <- c()
  i3_all <- c();  j3_all <- c()
  
  for (k in 1:sim.N){
    fname <- paste("simresult_N_",numobs,"_",k,".Rdata",sep="")
    fname1 <- paste("otherresult_N_",numobs,"_",k,".Rdata",sep="")
    fnamep <- paste("paradata_N_",numobs,"_",k,".Rdata",sep="")
    fname.data <- paste("simdata_N_",numobs,"_",k,".Rdata",sep="")
    
    
    load(file=fname.data)#data2
    load(file=fname) #sim.result
    load(file=fname1) #other.result
    load(file=fnamep) #parameters
    
    sim.mucontrol <- paradata[,2]
    sim.pcontrol <- paradata[,3]
    sim.mucase <- paradata[,5]
    sim.pcase <- paradata[,6]
    sim.geneid <- paradata[,7]
    
    datatmp <- data2
    count <- rowSums(datatmp==0)
    rawpall  <- count/2*numobs  
    
    est <- sim.result$estimates
    
    if (testmethod=="Both"){
      p <-  sim.result$pvalue 
      t.p <- other.result[,1] 
      twot.p <- other.result[,24]
    }
    
    if (testmethod=="Mean"){
      p <-  sim.result$pvalue_mean 
      t.p <- other.result[,7]
      twot.p <- other.result[,23]
    }
    
    if (testmethod=="Pzero"){
      p <-  sim.result$pvalue_zero
      t.p <- other.result[,13]
      twot.p <- other.result[,22]
    }
    
    aft.p <- other.result[,19]
    
    q<-p.adjust(p,method='BH')
    t.q<-p.adjust(t.p,method='BH')
    aft.q <- p.adjust(aft.p,method="BH")
    twot.q <- p.adjust(twot.p,method="BH")
    
    rank <- c(1:length(sim.geneid))
    
    a1_ <- sum(q[!is.na(q)]<=0.01 & sim.geneid[!is.na(q)]==1)
    b1_ <- sum(q[!is.na(q)]<=0.01 & sim.geneid[!is.na(q)]==0)
    c1_ <- sum(t.q[!is.na(t.q)]<=0.01 & sim.geneid[!is.na(t.q)]==1)
    d1_ <- sum(t.q[!is.na(t.q)]<=0.01 & sim.geneid[!is.na(t.q)]==0)
    e1_ <- sum(aft.q[!is.na(aft.q)] <= 0.01 & sim.geneid[!is.na(aft.q)]==1)
    f1_ <- sum(aft.q[!is.na(aft.q)] <= 0.01 & sim.geneid[!is.na(aft.q)]==0)
    i1_ <- sum(twot.q[!is.na(twot.q)] <= 0.01 & sim.geneid[!is.na(twot.q)]==1)
    j1_ <- sum(twot.q[!is.na(twot.q)] <= 0.01 & sim.geneid[!is.na(twot.q)]==0)
    
    a2_ <- sum(q[!is.na(q)]<=0.05 & sim.geneid[!is.na(q)]==1)
    b2_ <- sum(q[!is.na(q)]<=0.05 & sim.geneid[!is.na(q)]==0)
    c2_ <- sum(t.q[!is.na(t.q)]<=0.05 & sim.geneid[!is.na(t.q)]==1)
    d2_ <- sum(t.q[!is.na(t.q)]<=0.05 & sim.geneid[!is.na(t.q)]==0)
    e2_ <- sum(aft.q[!is.na(aft.q)] <= 0.05 & sim.geneid[!is.na(aft.q)]==1)
    f2_ <- sum(aft.q[!is.na(aft.q)] <= 0.05 & sim.geneid[!is.na(aft.q)]==0)
    i2_ <- sum(twot.q[!is.na(twot.q)] <= 0.05 & sim.geneid[!is.na(twot.q)]==1)
    j2_ <- sum(twot.q[!is.na(twot.q)] <= 0.05 & sim.geneid[!is.na(twot.q)]==0)
    
    a3_ <- sum(q[!is.na(q)]<=0.1 & sim.geneid[!is.na(q)]==1)
    b3_ <- sum(q[!is.na(q)]<=0.1 & sim.geneid[!is.na(q)]==0)
    c3_ <- sum(t.q[!is.na(t.q)]<=0.1 & sim.geneid[!is.na(t.q)]==1)
    d3_ <- sum(t.q[!is.na(t.q)]<=0.1 & sim.geneid[!is.na(t.q)]==0)
    e3_ <- sum(aft.q[!is.na(aft.q)] <= 0.1 & sim.geneid[!is.na(aft.q)]==1)
    f3_ <- sum(aft.q[!is.na(aft.q)] <= 0.1 & sim.geneid[!is.na(aft.q)]==0)
    i3_ <- sum(twot.q[!is.na(twot.q)] <= 0.1 & sim.geneid[!is.na(twot.q)]==1)
    j3_ <- sum(twot.q[!is.na(twot.q)] <= 0.1 & sim.geneid[!is.na(twot.q)]==0)
    
    a1_all <- c(a1_all,a1_)
    b1_all <- c(b1_all,b1_)
    c1_all <- c(c1_all,c1_)
    d1_all <- c(d1_all,d1_)
    e1_all <- c(e1_all,e1_)
    f1_all <- c(f1_all,f1_)
    i1_all <- c(i1_all,i1_)
    j1_all <- c(j1_all,j1_)
    
    a2_all <- c(a2_all,a2_)
    b2_all <- c(b2_all,b2_)
    c2_all <- c(c2_all,c2_)
    d2_all <- c(d2_all,d2_)
    e2_all <- c(e2_all,e2_)
    f2_all <- c(f2_all,f2_)
    i2_all <- c(i2_all,i2_)
    j2_all <- c(j2_all,j2_)
    
    a3_all <- c(a3_all,a3_)
    b3_all <- c(b3_all,b3_)
    c3_all <- c(c3_all,c3_)
    d3_all <- c(d3_all,d3_)
    e3_all <- c(e3_all,e3_)
    f3_all <- c(f3_all,f3_)
    i3_all <- c(i3_all,i3_)
    j3_all <- c(j3_all,j3_)
    
    sim.sort.result <- cbind(sim.geneid,q,p)
    sim.sort.result <-sim.sort.result[order(p),]
    sim.sort.result <- cbind(sim.sort.result,rank)
    
    t.sort.result <- cbind(sim.geneid,t.q,t.p)
    t.sort.result <-t.sort.result[order(t.p),]
    t.sort.result <- cbind(t.sort.result,rank)
    
    t2.sort.result <- cbind(sim.geneid,aft.q,aft.p)
    t2.sort.result <-t2.sort.result[order(aft.p),]
    t2.sort.result <- cbind(t2.sort.result,rank)
    
    twot.sort.result <- cbind(sim.geneid,twot.q,twot.p)
    twot.sort.result <-twot.sort.result[order(twot.p),]
    twot.sort.result <- cbind(twot.sort.result,rank)
    
    #tOP RANKED FEATURES AND TRUE POSITIVE RATE   
    TOPn <- seq(5,TOPNMAX,5)
    TPR <- c()
    for (i in TOPn){
      resultsub <- sim.sort.result[1:i,]
      DE <- resultsub[,1][resultsub[,1]==1]
      DE[is.na(DE)] <- 0
      tpr <- length(DE)/nrow(resultsub)  
      TPR <- c(TPR,tpr)
      TPR[is.na(TPR)] <- 0
    }
    DASEV.TPR <-cbind(DASEV.TPR,TPR)
    
    
    T.TPR <- c()
    for (i in TOPn){
      resultsub <- t.sort.result[1:i,]
      DE <- resultsub[,1][resultsub[,1]==1]
      DE[is.na(DE)] <- 0
      tpr <- length(DE)/nrow(resultsub)  
      T.TPR <- c(T.TPR,tpr)
      T.TPR[is.na(T.TPR)] <- 0
    }
    TLK.TPR <- cbind(TLK.TPR,T.TPR)
    
    
    T2.TPR <- c()
    for (i in TOPn){
      resultsub <- t2.sort.result[1:i,]
      DE <- resultsub[,1][resultsub[,1]==1]
      DE[is.na(DE)] <- 0
      tpr <- length(DE)/nrow(resultsub)  
      T2.TPR <- c(T2.TPR,tpr)
      T2.TPR[is.na(T2.TPR)] <- 0
    }
    AFT.TPR <- cbind(AFT.TPR,T2.TPR)
    
    
    TWOT.tpr <- c()
    for (i in TOPn){
      resultsub <- twot.sort.result[1:i,]
      DE <- resultsub[,1][resultsub[,1]==1]
      DE[is.na(DE)] <- 0
      tpr <- length(DE)/nrow(resultsub)  
      TWOT.tpr <- c(TWOT.tpr,tpr)
      TWOT.tpr[is.na(TWOT.tpr)] <- 0
    }
    TWOT.TPR <- cbind(TWOT.TPR,TWOT.tpr)
    
}
  
  mean.DASEV.TPR<-rowMeans(DASEV.TPR)
  mean.TLK.TPR <-rowMeans(TLK.TPR)
  mean.AFT.TPR<-rowMeans(AFT.TPR)
  mean.TWOT.TPR <- rowMeans(TWOT.TPR)
  
  a1<- round(mean(a1_all))
  b1<- round(mean(b1_all))
  c1<- round(mean(c1_all))
  d1<- round(mean(d1_all))
  e1<- round(mean(e1_all))
  f1<- round(mean(f1_all))
  i1<- round(mean(i1_all))
  j1<- round(mean(j1_all))
  
  
  a2<- round(mean(a2_all))
  b2<- round(mean(b2_all))
  c2<- round(mean(c2_all))
  d2<- round(mean(d2_all))
  e2<- round(mean(e2_all))
  f2<- round(mean(f2_all))
  i2<- round(mean(i2_all))
  j2<- round(mean(j2_all))
  
  a3<- round(mean(a3_all))
  b3<- round(mean(b3_all))
  c3<- round(mean(c3_all))
  d3<- round(mean(d3_all))
  e3<- round(mean(e3_all))
  f3<- round(mean(f3_all))
  i3<- round(mean(i3_all))
  j3<- round(mean(j3_all))
  
  
  v1 <- c(a1,a2,a3)
  v2 <- c(-b1,-b2,-b3)
  v3 <- c(c1,c2,c3)
  v4 <- c(-d1,-d2,-d3)
  v5 <- c(e1,e2,e3)
  v6 <- c(-f1,-f2,-f3)
  v9 <- c(i1,i2,i3)
  v10 <- c(-j1,-j2,-j3)
  tfdrd1 <- b1_all/(b1_all+a1_all)
  tfdrd2 <- b2_all/(b2_all+a2_all)
  tfdrd3 <- b3_all/(b3_all+a3_all)
  tfdrt1 <- d1_all/(d1_all+c1_all)
  tfdrt2 <- d2_all/(d2_all+c2_all)
  tfdrt3 <- d3_all/(d3_all+c3_all)
  tfdra1 <- f1_all/(f1_all+e1_all)
  tfdra2 <- f2_all/(f2_all+e2_all)
  tfdra3 <- f3_all/(f3_all+e3_all)
  tfdrtwot1 <- j1_all/(j1_all+i1_all)
  tfdrtwot2 <- j2_all/(j2_all+i2_all)
  tfdrtwot3 <- j3_all/(j3_all+i3_all)
  
  tprrate <- c(mean(tfdrd1[!is.na(tfdrd1)]),mean(tfdrt1[!is.na(tfdrt1)]),
               mean(tfdra1[!is.na(tfdra1)]),mean(tfdrtwot1[!is.na(tfdrtwot1)]),
               mean(tfdrd2[!is.na(tfdrd2)]),mean(tfdrt2[!is.na(tfdrt2)]),
               mean(tfdra2[!is.na(tfdra2)]),mean(tfdrtwot2[!is.na(tfdrtwot2)]),
               mean(tfdrd3[!is.na(tfdrd3)]),mean(tfdrt3[!is.na(tfdrt3)]),
               mean(tfdra3[!is.na(tfdra3)]),mean(tfdrtwot3[!is.na(tfdrtwot3)]))
  rate <- paste(round(100*tprrate, 2), "%", sep="")
 
  res <- list(TOPn=TOPn,DASEV.TPR=mean.DASEV.TPR,TLK.TPR=mean.TLK.TPR,
              AFT.TPR= mean.AFT.TPR,TWOT.TPR=mean.TWOT.TPR,
              v1=v1,v2=v2,v3=v3,v4=v4,v5=v5,v6=v6,v9=v9,v10=v10,rate=rate)

}
######################################################################

##########################SECOND PART#################################
#################Call Function to get average results#################
######################################################################
TOPNMAX =300
##########################Secnario 1##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S1/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S1/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S1/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S1/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S1N10.Rdata")
save(res2,file="S1N20.Rdata")
save(res3,file="S1N100.Rdata")
save(res4,file="S1N200.Rdata")
######################################################################

##########################Secnario 2##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S2/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S2/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S2/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S2/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S2N10.Rdata")
save(res2,file="S2N20.Rdata")
save(res3,file="S2N100.Rdata")
save(res4,file="S2N200.Rdata")
######################################################################

##########################Secnario 3##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S3/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S3/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S3/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S3/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S3N10.Rdata")
save(res2,file="S3N20.Rdata")
save(res3,file="S3N100.Rdata")
save(res4,file="S3N200.Rdata")
######################################################################

##########################Secnario 4##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S4/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S4/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S4/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S4/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S4N10.Rdata")
save(res2,file="S4N20.Rdata")
save(res3,file="S4N100.Rdata")
save(res4,file="S4N200.Rdata")
######################################################################

##########################Secnario 5##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S5/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S5/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S5/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S5/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S5N10.Rdata")
save(res2,file="S5N20.Rdata")
save(res3,file="S5N100.Rdata")
save(res4,file="S5N200.Rdata")
######################################################################

##########################Secnario 6##################################
######################################################################
res1 <- Figure4f(dirpath ="c:/DASEV/S6/", testmethod = "Both",
                 sim.N =100,numobs=10,TOPNMAX=TOPNMAX)

res2 <- Figure4f(dirpath ="c:/DASEV/S6/", testmethod = "Both",
                 sim.N =100,numobs=20,TOPNMAX=TOPNMAX)

res3 <- Figure4f(dirpath ="c:/DASEV/S6/", testmethod = "Both",
                 sim.N =100,numobs=100,TOPNMAX=TOPNMAX)

res4 <- Figure4f(dirpath ="c:/DASEV/S6/", testmethod = "Both",
                 sim.N =100,numobs=200,TOPNMAX=TOPNMAX)
setwd("c:/DASEV/Figure4/")
save(res1,file="S6N10.Rdata")
save(res2,file="S6N20.Rdata")
save(res3,file="S6N100.Rdata")
save(res4,file="S6N200.Rdata")
######################################################################

###########################THIRD PART#################################
#############################Plots####################################
######################################################################
setwd("c:/DASEV/Figure4/")
xlimmax=300
setcex <- 1
par(mfrow=c(2,4))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))

############################Figure S8#################################
load("S1N10.Rdata") #N=10 res1
load("S1N20.Rdata") #N=20 res2
load("S1N100.Rdata") #N=100 res3
load("S1N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)-70), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)), labels=res$rate,cex=setcex)
######################################################################

############################Figure S9#################################
load("S2N10.Rdata") #N=10 res1
load("S2N20.Rdata") #N=20 res2
load("S2N100.Rdata") #N=100 res3
load("S2N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)-70), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)), labels=res$rate,cex=setcex)
######################################################################

############################Figure S10################################
load("S3N10.Rdata") #N=10 res1
load("S3N20.Rdata") #N=20 res2
load("S3N100.Rdata") #N=100 res3
load("S3N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,210,250,210,250,250,290,250,290)-70), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,180,150,180,200,230,200,230,220,245,220,245))+50, labels=res$rate,cex=setcex)
######################################################################

############################Figure S11################################
load("S4N10.Rdata") #N=10 res1
load("S4N20.Rdata") #N=20 res2
load("S4N100.Rdata") #N=100 res3
load("S4N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,210,250,210,250,250,290,250,290)-70), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,180,150,180,200,230,200,230,220,245,220,245))+50, labels=res$rate,cex=setcex)
######################################################################

############################Figure S12################################
load("S5N10.Rdata") #N=10 res1
load("S5N20.Rdata") #N=20 res2
load("S5N100.Rdata") #N=100 res3
load("S5N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,400),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,210,250,210,250,250,290,250,290)-30), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,400),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,180,150,180,200,230,200,230,220,245,220,245))+130, labels=res$rate,cex=setcex)
######################################################################

############################Figure S13################################
load("S6N10.Rdata") #N=10 res1
load("S6N20.Rdata") #N=20 res2
load("S6N100.Rdata") #N=100 res3
load("S6N200.Rdata") #N=200 res4

res <- res1
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
legend("topright", legend=c("DASEV","TLK","AFT","2T"),border="white", col=c("darkorange2","royalblue4","darkred","seagreen"),lty=c(1,1),cex=1.5)
res <- res2
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res3
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="c",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res4
plot(res$TOPn,res$DASEV.TPR ,type="l",col="darkorange2", xlim=c(0,xlimmax),ylim=c(0,1), xlab="Top ranked features", ylab="True positive rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="d",lwd=2)
lines(res$TOPn,res$TLK.TPR,type="l",col="royalblue4", lwd=2)
lines(res$TOPn,res$AFT.TPR,type="l",col="darkred", lwd=2)
lines(res$TOPn,res$TWOT.TPR,type="l",col="seagreen", lwd=2)
res <- res1
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "e",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
legend("bottomleft",ncol=2, legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP","AFT TP","AFT FP","2T TP","2T FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2","darkred","pink","seagreen","lightseagreen"),cex=1.5,pch=c(15,15,15,15,15,15,15,15))
res <- res2
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "f",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,150,190,150,190,150,190,150,190)-90), labels=res$rate,cex=setcex)
res <- res3
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "g",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)-70), labels=res$rate,cex=setcex)
res <- res4
res$rate[res$rate=="NaN%"] <- "N.A."
barplot(matrix(c(res$v1,res$v3,res$v5,res$v9),nrow=4, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4',"darkred","seagreen"),ylim=c(-200,300),xlab="Reported FDR", ylab = "Counts",main = "h",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(res$v2,res$v4,res$v6,res$v10),nrow=4, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2',"pink","lightseagreen"), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,3.5,4.5,6.5,7.5,8.5,9.5,11.5,12.5,13.5,14.5), y=(c(150,190,150,190,180,220,180,220,210,250,210,250)), labels=res$rate,cex=setcex)
######################################################################